Optimal Power Flow

Unlike the ED model, the optimal power flow (OPF) takes into account network constraints. The latter can be achieved by using either a dc-based or ac-based power flow model. In the following, we formulate and solve both the dc-based OPF model for a simple 3 bus example using Julia's JuMP package.

The objective function of the OPF model can be formulated as follows: $$ \min \sum_{i \in I} c^g_{i} \cdot g_{i} + \sum_{b \in B} c^w \cdot w_b, $$ where $b$ is the index of buses. Note that as compared to the ED model, the OPF accounts for wind power injection at each bus separately.

s.t.

Minimum ($g^{\min}$) and maximum ($g^{\max}$) limits on power outputs of generators:

$$ g^{\min}_{i} \leq g_{i} \leq g^{\max}_{i}. $$

Constraint on the wind power injection:

$$ 0 \leq w_b \leq w^f_b, \quad \forall b \in B $$ where $w_b$ and $w^f_b$ are the nodal wind power injection and wind power forecast, respectively.

Nodal power balance constraint:

$$ \sum_{i \in I_b} g_{i} + w_b + \sum_{l \in L_b} f_l= d_b^f, \quad \forall b \in B $$ where $d^f_b$ is the nodal demand forecast, $f_l$ is the power flow in line $l$, and $L_b$ is the set of lines connected to bus $b$.

DC power flow model

Power flow limit on every transmission line:

$$ -f^{\max}_{l} \leq f_{l} = \frac{\Delta \theta_{l}}{x_{l}} \leq f^{\max}_{l}, \quad \forall l \in L, $$ where $f_{l}$ and $f^{\max}_{l}$ are the power flow and the maximum power flow limit of line $l$, respectively; and $x_{l}$ and $\Delta \theta_{l}$ are the angular difference between the ends of line $l$ and its reactance, respectively.

Angular difference:

$$ -\pi \leq \Delta \theta_{l} \leq \pi, \quad \forall l \in L. $$

Implementation of DC Optimal Power Flow



In [ ]:

    
using JuMP
using Interact
using Gadfly
using GraphViz

In the following experiments we employ the three bus example as in the ED and UC experiments.



In [ ]:

    
# Maximum power output of generators (card(I))
g_max = [1000 1500];
# Minimum power output of generators (card(I))
g_min = [0 300];
# Incremental cost of generators (card(I))
c_g = [50 100];
# Generator map (card(I) x card(b))
g_map = [0 1 0; 0 0 1];
# Incremental cost of wind generators
c_w = 50;
# Total demand (card(B))
d = [0 0 1500];
# Wind forecast (card(V))
w_f = [150 50];
# Generator map (card(I) x card(b))
w_map = [1 0 0; 0 1 0];
# Power flow limits (card(L))
f_max = [100 1000];
# Line map (card(L) x card(B))
f_map = [1 -1 0; 0 1 -1]; 
# Line impedance 
x = [0.001 0.001]

Sensitivity of the OPF solution to power flow limits

To assess the impact of power flow limits, we will use the manipulator, varying the power flow in line L2, and observe the distribution of power flows among transmission lines. This graphical representation is enabled by package GraphViz.



In [ ]:

    
function writeDot(name, busIdx, busInj, renGen, f, t, lineFlow, lineLim, size=(11,14))
    # This function generates a graph that richly expresses the RTS96 system state.
    # name              a name for the graph and resulting dot file
    # busIdx            bus names (could be text) in order of busInj
    # busInj            injection at each bus
    # renGen            renewable generation at each bus (0 for non-wind buses)
    # f                 "from" node for each line
    # t                 "to" node for each line
    # lineFlow          flow on each line
    # lineLim           limit for each line
    # size              size of graphical output

    busInj = round(busInj,2)
    lineFlow = round(lineFlow,2)
    
    # Open the dot file, overwriting anything already there:
    dotfile = IOBuffer()
    
    # Begin writing the dot file:
    write(dotfile, "digraph $(name) {\nrankdir=LR;\n")

    # Set graph properties:
    write(dotfile, 
    "graph [fontname=helvetica, tooltip=\" \", overlap=false, size=\"$(size[1]),$(size[2])\", ratio=fill, orientation=\"portrait\",layout=dot];\n")

    # Set default node properties:
    write(dotfile, "node [fontname=helvetica, shape=square, style=filled, fontsize=20, color=\"#bdbdbd\"];\n")

    # Set default edge properties:
    write(dotfile, "edge [fontname=helvetica, style=\"setlinewidth(5)\"];\n")

    # Write bus data to dot file:
    for i = 1:length(busIdx)
        write(dotfile, 
        "$(i) [label=$(int(busIdx[i])), tooltip=\"Inj = $(busInj[i])\"") # bus label and tooltip

        # Represent renewable nodes with blue circles:
        if union(find(renGen),i) == find(renGen)
            write(dotfile, ", shape=circle, color=\"#5677fc\"")
        end

        write(dotfile, "];\n")
    end

    # Write line data to file:

    for i = 1:length(f)

        normStrain = abs(lineFlow[i])/lineLim[i] # normalized strain on line i

        # Use flow direction to determine arrow direction,
        # label with flow,
        # and color according to strain
        if lineFlow[i] > 0
            write(dotfile, 
            "$(f[i]) -> $(t[i]) [label=$(lineFlow[i])")
        else
            write(dotfile, 
            "$(t[i]) -> $(f[i]) [label=$(-lineFlow[i])")
        end
        write(dotfile,
        ", tooltip=\" \", labeltooltip=\"Flow = $(int(normStrain*100))%\", color=\"$(round((1 - normStrain)/3,3)) 1.000 0.700\"];\n")
    end

    write(dotfile, "}\n")

    dottext = takebuf_string(dotfile)
    #print(dottext)

    return Graph(dottext)
end



In [ ]:

    
@manipulate for line2_limit = 1000:10:1500
    f_max[2] = line2_limit
#function solve_dcopf (g_max, g_min, c_g, c_w, d, w_f)
#Define the optimal power flow (OPF) model
opf=Model() 
    
# Define decision variables    
@defVar(opf,   g[i=1:2] >= 0 ) ; # power output of generators
@defVar(opf, w[v=1:2] >=0 ) ; # wind power injection
@defVar(opf, f[l=1:2]) ; # power flows 
@defVar(opf, theta[b=1:3]) ; # bus angle 


# Define the objective function
@setObjective(opf,Min,sum{c_g[i] * g[i],i=1:2} + sum{c_w * w[v],v=1:2});


for i in 1:2
    @addConstraint(opf,  g[i] <= g_max[i] ) #maximum
    @addConstraint(opf,  g[i] >= g_min[i] ) #minimum
end


# Define the constraint on the wind power injections
for v in 1:2
    @addConstraint(opf,  w[v] <= w_f[v]); 
end

# Define the constraint on the power flows
for l in 1:2
    @addConstraint(opf, f[l] <= f_max[l]); # direct flows
    @addConstraint(opf, f[l] >=  -f_max[l]); # reverse flows
end

# Define the power balance constraint
for b in 1:3
    @addConstraint(opf, sum{g_map[i,b]* g[i],i=1:2} + sum{w_map[v,b] * w[v], v=1:2} + sum{f_map[l,b] * f[l], l=1:2}>= d[b]); 
end

# Calculate f[l]
for l in 1:2
    @addConstraint(opf, f[l] == 1/x[l] * sum{f_map [l,b] * theta[b], b=1:3})  ; # power flow in every line
end

# Slack bus 
    @addConstraint(opf, theta [b=1] == 0)  ; # direct flows


# Solve statement
solve(opf)
writeDot("UC",[1,2,3],[0.0,getValue(g)[:]],[getValue(w)[:],0.0],[2,3],[1,2],getValue(f)[:],f_max,(3,3))
end



In [ ]: